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ABSTRACT 



Double-diffusive convection, often referred to as semi-convection in astro- 
pliysics, occurs in thermally and compositionally stratified systems which are sta- 
ble according to the Ledoux-criterion but unstable according to the Schwarzchild 
criterion. This process has been given relatively little attention so far, and its 
properties remain poorly constrained. In this paper, we present and analyze a set 
of three-dimensional simulations of this phenomenon in a Cartesian domain under 
the Boussinesq approximation. We find that in some cases the double-diffusive 
convection saturates into a state of homogeneous turbulence, but with turbulent 
fluxes several orders of magnitude smaller than those expected from direct over- 
turning convection. In other cases the system rapidly and spontaneously develops 
closely-packed thermo-compositional layers, which later successively merge until 
a single layer is left. We compare the output of our simulations with an existing 
theory of layer formation in the oceanographic context, and find very good agree- 
ment between the model and our results. The thermal and compositional mixing 
rates increase significantly during layer formation, and increase even further with 
each merger. We find that the heat fiux through the staircase is a simple func- 
tion of the layer height. We conclude by proposing a new approach to studying 
transport by double-diffusive convection in astrophysics. 



Subject headings: convection - hydrodynamics - planets and satellites: general - 
stars:interior 
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Introduction 



1.1. Convection, double-diffusive convection (semi-convection) and fingering 

convection (thermohafine convection) 

One of tlie longest- standing problems in stellar and planetary astrophysics is that of 
modeling the transport of heat and chemical species within turbulent regions. The best- 
studied and most ubiquitously relevant case is that of overturning convection through a 
chemically homogeneous gas layer. There, the well-known Schwarzchild criterion is used to 
determine the extent of the convective region, while the transport prope rties through the 
layer are commonly modeled using mixing-length theory ( iBiermannI 119321 ) . The success of 
well-calibrated mixing-length models in explaining many observable properties of stars is 
quite remarkable. 

However, much less is known about convection in the presence of addition al factors suc h 



as strong rotation, strong magnetic fields and strong compositional g radients (ISpiege 



In all cases, the linear stability of the system is well-understood (j Chandr asekhar 



1972 ). 



196l|) 



but characterizing its fully-nonlinear transport properties remains the subject of ongoing 
research. In this work, we focus on the case of convection in the presence of a strong stabi- 
lizing compositional gradient, but in the absence of rotation or magnetic field . This regime 
is often called "semi-convection" in astrophysics (jSchwarzschild fc Harmlll958l ). although we 
prefer to use the terminology "double-diffusive convection" commonly used in oceanography 
to clarify the true nature of the instability responsible for the turbulence. 

It has long been known that the relevant criterion for instability to overturning convec- 
tion in the presence of a compositional gradient is not the Schwarzchild criterion. 
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where T is the temperature, /i the mean molecular weight, p the pressure, and where the 
subscript "ad" expresses a derivative at constant specific entropy. In fact, both of these 
criteria merely express the same property when written in terms of the density stratification: 



(3) 
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A system is unstable to overturning convection if the density of a parcel of fluid, raised 
adiabatically and in pressure equilibrium from its original position, is lower than that of its 
new surroundings. 

The question of what happens to regions which are stable according t o the Ledoux crite- 
rion b ut unstable according to the Schwarzchild criterion was first raised bv lSchwarzschild fc Harm 



(119581) . It was later found that this regime is in fact also linearly unstable (jWalirull964j : iKato 
19661 ). but through a double- diffusive instability, i.e. an instability which cannot occur un- 
less the thermal diffusivity of the fluid, k^, is larger than its compositional diffusivity k^. 
This condition is however automatically satisfied in stellar and planetary interiors where the 
diffusivity ratio r = k^/kt (often called the inverse Lewis number), can be as low as 10~^. 

As a result, a wide range of situations arise in which double-diffusive convection occurs 
and controls transport within the ob ject. A commonly studied case is that of semi-convection 



at the edge of co re-convective stars (lLedouxlll947l : lTayleiill954l : ISchwarzschild &: Harmlll958 



Merryfieldlll995l ). In moderately massive stars for example, a mean molecular weight gradi- 
ent develops over time at the edge of the core, and eventually begins to affect convection. 
When it is strong enough to stabilize the fluid, the fully-convective region shrinks in size, 
leaving behind a "semi-convective" region in which transpor t is controlled by double- diffusive 
processes instead. Other related examples are reviewed by iMerryfieldl ( 119951 ). 



The possible role of double-diffusive convection in regulating thermal and composi- 
tional transport in the interior of giant planets was recognized later, and has been discussed 
i n the context of convective planetary envelopes where the stabilizing component is helium 
(IStevenson fc Salpeterl 119771 ) . higher- metallicity material at the edge of a rocky core (see 
StevensonI (119851 ) and in particular Figure 2 of his paper), or methane (IGierasch fc Conrath 
19871 ). Double-diffusive convection ha s also been proposed to ex plain the abnormally large 
radii of some transiting exoplanets ( IChabrier fc Baraffd 120071 ). Finall y, it has recently 
been invoked as a new mechanism for driving pulsations in white dwarfs (jShibahashil 120071 : 
Kurtz et al.ll2008l ). 



Before we move on to describe existing work on double-diffusive convection, we note 
that it should n ot be confus ed with that arising from the related double-diffusive "finger- 
ing" instability ( 1Sternlll960l ). The latter also occurs in Ledoux-stable systems, but in the 
opposite situation when the more rapidly diffusing thermal field is stably stratified while 
the slowly- diffusing compositional field is unstably stratified. Its turbulent manifestation is 
often referred to as "thermohaline convection" in astrophysics, by analogy with the oceanic 
case in which the compositional gradient is due to salt. We prefer a more general terminol- 
ogy and use the alternative name "fingering c onvection". F i gure [2 illustrates for clarity the 
various regimes of convective instability. See iTraxler et al.l ( l2010al ) for a study of fingering 
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convection in the astrophysical regime. 



1.2. Previous work on double-diffusive convection 



Very little is known about mixing by double-diffusive convection, despite its obvious 
importance in stellar and planetary astrophysic s. Linear stability reveals that the unstabl e 
modes take the form of overstable gravity waves JWaliJl964l : lKatJl966l : ISaines fc Gillll969h . 
What governs the saturation of the instability in the astrophysically-relevant parameter 
regime, however, remains essentially unknown. Note that by "astrophysically-relevant" we 
imply a low diffusivity ratio, r -C 1, and a low Prandtl number, Pr = z//kt ^ 1, where v is 
the viscosity. Both numbers are typically of the order of 10~^ — 10~^ in stellar and planetary 
interiors. 

To add to the complexity of the problem, double- diffusive convection is known in some 
cases to lead to thermo-compositional layering, i.e. to the development of stacks of well-mixed 
fuUy-convective layers separated by strongly stratified int erfaces. This double-diffusive layer- 
ing is commonly observed for example in the arctic ocean (INeal et al.lll969l : iToole et al.ll2006[ ) 
where cool fresh water lies on top of warmer, saltier water. It has been studied exteiisively 



in laboratory experiments ( lTurnerlll968l : iHuppert &: Lindenlll979l : iHuppert fc Turner! Il980 



Turner! Il985l ). An important result of these studies is that turbulent mixing in the presence 
of layers is significantly enhanced compared with that of a system which has the same overall 
contrast in temperature and composition, but where the stratification is everywhere much 
smoother. 

Whether layer formation occurs in double-diffusive co nvection at l ow Prandtl number 
actually remains to be determined - it is usually assumed ( lSpruit!ll992l : IChabrier fc Baraffe 
20071 ). by analogy with the high-Prandtl number oceanographic case. It is important to 
realize, howev e r, that such analogies can be misleading. This was recently demonstrated by 
Traxler et al.l ( l2010a( ) in the case of fingering convection. Similar th ermohaline staircase s 
are ubiquitously observed in fingering-unstable regions of the ocean ( ISchmitt et al.!l2005l ). 
and have been shown to form spontaneously through a secondary linear instability of h omo- 
geneous fingering convection (IRadko! l2003l : iTraxler et al.! 12010b! : IStellmach et al.! 120101 ) (see 
^too). However, iTraxler et al.! ( l2010al ) demonstrated that this secondary instability does 
not happen in the astrophysical context and concluded that thermo-compositional layers are 
not expected in that case. In other words, given that the analogy with the heat-salt system 
doesn't hold in the fingering regime, one should be extremely cautious about using it a priori 
in the double-diffusive regime. 
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To summarize, it is known that "homogeneous" and "layered" double-diffusive convec- 
tion have rather different mixing properties. A good mixing parametrization therefore needs 
to incorporate both cases, and must include a criterion to decide whether the system con- 
sidered lies in one or the other. Exist ing parametrizations, by con t rast, have so far ei ther 
ignored the possible effe ct of layering (jSchwarzschild fc HarmI Il958l : iLanger et al.l 119831 ) , or 
relied on it (ISpruitI 119921 ). 



A few numerical simulations of double-diffusive convection have been performed to 
date to address t he problem. The first of this kind (to our knowledge), were presented by 
Merryfieldl (Il995[ ). He ran a series of two-dimensional (2D) anelastic simulations, horizontally- 
periodic, and bounded in the vertical direction by two plates. Double-diffusive convection 
was forced through the imposed boundary conditions, whi ch maintained a g iven overall tem- 
perature and compositional contrast across the domain. iMerryfieldl (119951 ) focused on un- 
derstanding how the mechanism responsible for the saturation of the initial double-diffusive 
instability depends on the governing parameters, and in particular the strength of the ther- 
mal driving, the Prandtl number and the diffusivity ratio. He also compared the outcome 
of his simulations with existing parametrizations of double- diffusive convection both in the 
absence of layers, and in the presence of initially forced layers. One of the main difficul- 
ties encountered was the development of numerical instabilities in cases with strong thermal 
driving, which prevented him from drawing definite conclusions on the long-term statistical 
properties of the turbulence. In addition, in many of the runs the flows were eventually 
influenced by the presence of domain boundaries. 

In sub sequent years , two additional attempts at modeling double-diffusive convection 
were made. iBielld ( 1200 ll ). as part of his PhD thesis, ra n a series of 2D fully- compressible sim- 
ulations which complement those of IMerryfieldl ( 1l995l ). There, the system was also confined 
between two plates, but the boundary conditions were "fixed fiux" conditions. iBielld (120011 ) 
was interested in studying more specifically the layer formation process, and his experimen 



tal se tup was similar to that of heat-salt laboratory experiments (e.g. iHuppert fc Linden 
19791 ). For this purpose, the simulations were initialized with stable uniform gradients in 
temperature and composition, but destabilized at t = by increasing the heat fiux at the 
bottom boundary. He found that the first bottom layer easily forms, but did not observe any 
subsequent layer formation. He analyzed the dynamics of the interface, and concluded that 
interfacial transport was dominated by wave-breaking rather than by diffusive processes as 
is often assumed. Howeve r, boundary eff ects also began to influence the results of his sim- 
ulations after some time. iBascoull ( 120071 ). also as part of his PhD thesis, studied a similar 
2D time- dependent system, in which the initial background state had a homogeneous com- 
position and neutrally stable temperature gradient, and where the system was destabilized 
by an imposed heat and mean-molecular weight flux through the bottom boundary. He 
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also observed the growth of a convective layer near the bottom boundary and studied its 
development, for high Prandtl number (heat-salt regime) and for low Prandtl number ("as- 
trophysical" regime). He was unable, however, to run his simulations long enough to achieve 
statistical equilibrium. 

In this paper, we present a new series of three-dimensional (3D) numerical experiments 
to study mixing by double-diffusive convection. We approach the problem from a different 
but complementary angle, and try to address some of the inherent shortcomings of the 
experimental setup used in previous studies: we focus specifically on measuring the quasi- 
steady statistical properties of double- diffusive turbulence, and use a numerical setup which 
minimizes the effects of domain boundaries. 

We discuss the model setup and briefly summarize its linear stability properties in §2J 
The numerical algorithm and the selection of the experimental parameters are described 
in ^ In each described in §H we extract values for the transport coefficients 

while the system is in a state of homogeneous turbulence. However, we find that for more 
unstable systems a secondary instability leads to the formation of thermo-compositional 
layers. The layers continue evolving and successively merge, and each merger is accompanied 
with a significant increase in the transport coefficients. These results are discussed in_§5l 



and compared with a recent theory of layer formation in the oceanographic context (iRadko 



20031 ) . We find good agreement between this theory and our numerical results, which enables 
us to deduce a general criterion for the spontaneous formation of layers in double-diffusive 
convection in the astrophysical context. We discuss our results in §H1 



2. Model description and linear analysis 



In this section we present the governing equations and briefly summarize known results 
on the linear sta b ility of the problem. Our formalism is overall very similar to the one used 
by lTraxler et al.l (l2010al ) for fingering convection. 



2.1. Governing equations and boundary conditions 

As we demonstrate later (see §2.2p . the typical lengthscale of the unstable motions (in 
the absence of layers) is of the order of a few kilometers at most in parameter regimes 
typical of stellar and planetary interiors. It is therefore justified to neglect the effect of 
curvature entirely and work in a local Cartesian reference frame {x,y,z). Here, gravity 
defines the vertical direction: g = —gSz- We ignore the possible presence of magnetic fields 
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for simplicity, and neglect the effect of rotation. The latter is justified whenever the mean 
rotation rate is much smaller than the buoyancy (Brunt- Vaisala) frequency, which is often 
the case. 

In all that follows, we simplify the problem by using the Boussinesq approximation 



( ISpiegel fc Veronislll960l ). This approximation is a regular asymptotic limit of the primitive 
governing equations (for mass, momentum and energy conservation) when the domain height 
Lz is much smaller than the density scaleheight Dp (i.e. Lz/Dp — )■ 0), and when the typical 
flow velocity u is much smaller than the sound speed (i.e. ujcg — )■ 0). It is therefore 
particularly relevant here since the typical lengthscale of the perturbations is much smaller 
than a pressure or density scaleheight (unless the region considered is very close to the 
photosphere), and the typical velocities are always substantially subsonic (unless the system 
is very close to the onset of overturning convection). 

Since we consider a small fluid region within a much larger system, we can assume that 
the background composition and temperature gradients /zoz, Toz and Tq^ are constant. Here, 
the index denotes a background field, and the index z denotes a derivative with respect 
to the vertical coordinate z. We restrict our study to the case of double-diffusive convection 
by choosing the background stratification such that > V — Vad > 0, or equivalently 
/io. < To, - Tof < 0. 

In the Boussinesq approximation the mass conservation equation is replaced by the 
continuity equation 

V ■ u = , (4) 
where u = (w, f , w ) is the veloc ity field, and the (dimensional) thermal energy equation is 



approximated (e.g. lUlrichI 119721 ) by 



dT 

— + u ■ VT + (To, - T^)w = ktVT , (5) 

where T now represents the dimensional temperature perturbation (all background quanti- 
ties being denoted by the subscript instead), and we have assumed for simplicity that the 
thermal diffusivity is constant. This equation is derived from an energy conservation princi- 
ple, noting that in the Boussinesq approximation temperature and entropy perturbations are 
proportional. The term To, — therefore models the advection of the background entropy 
gradient. 

In the Boussinesq approximation, the density, temperature and mean molecular weight 
(p, T and respectively) are related via 

P- = -aT + , (6) 
Po 
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where po is the mean density of the region considered, and a and /3 are the coefficients of 
thermal expansion and compositional contraction respectively (e.g. a = 1/Tq and P = 1/fiQ 
for a perfect gas, where Tq and fiQ are the mean temperature and mean-molecular weight in 
the region considered). 

We construct our numerical model in such a way as to minimize the effects of the 
computational domain boundaries. For this purpose we use a triply-periodic box of size 
{L^, Ly, Lz), in which convection is permanently forced by the aforementioned background 
stratification. This approach has recently been used with success in mod eling and stud ying 



the formation of thermohaline staircases in the oceanographic context by iRadkd (|2003[ ) and 



Stellmach et al.l (|2010[ ). and is discussed in more detail in these papers. In this framework, 
the temperature and mean molecular weight fields can be written as the sum of a background 
variation plus perturbations which are triply-periodic functions of (x, y, z) such that 

T{x, y, z, t) = T{x + L^, y, z, t) = T{x, y + Ly, z, t) = T{x, y,z + L^, t) , (7) 

and similarly for fi. The pressure perturbation p and velocity field u are also assumed to be 
triply-periodic functions in the same way. 

We non-dimensionalize the equations using the anticipated lengthscale of the fastest 



grow ing modes of linear instability, which is a thermal diffusion scale (e.g. iBaines fc Gill 



19691): 

d 



ag\T,-,-T§t\) \m 

where N is the buoyancy frequency. Very roughly, in typical stellar interiors, u = O(10)cm^/ s, 
kt = O(10^)cm^/s and A^^ = O(10~^)s~^ so d = O(10^'^)cm, or in other words, a few hun- 
dreds of meters only. 

Note that with this definition, the thermal Rayleigh number defined on the finger scale 
is exactly one, while the global Rayleigh number. 



kti^ \ d 

is a function of the dimensionless height of the domain only. The unit timescale is taken to 
be the diffusion timescale across d, namely [t] = d'^/nx, and the velocity scale is [v] = Kx/d. 
The unit temperature is [T] = d\TQz — Tq^|, and the unit mean-molecular weight is [fi] = 
(a//3)|To2 — Tq^|c?. The resulting non-dimensional governing equations are 

1 / Ou \ 



Pr V dt 
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dT ~ 

— + U-VT -w = V^T , 
ot 

V-fi = 0, (10) 

where quantities with tildes now represent the dimensionless, triply-periodic perturbations. 
The non-dimensionalization introduces three parameters, namely the aforementioned Prandtl 
number Pr= u/kt and diffusivity ratio r = k^^/kt, as well as the so-called density ratio 

a|To,-Tof| V-Vad .... 

Finally, note that for reasons described in §2.31 it is common and preferable to work with 
the inverse density ratio Rq^ = I/Rq as a governing parameter. 



2.2. Linear stability 



Th e linear stability of this problem i s well- understood, thanks to the works of IWalin 



Veronisi (119651 ) and iBaines fc Gilll (119691 ) in the oceanographic context, and 



Kato 



19661 ) in the astrophysical context. The salient points are summarized here for completeness 



and clarity, as some of them will be used later. 

To analyze the stability of the governing equations we first linearize them around T 



ilx+imy+ikz+\t 



jl = and u = 0, and assume normal forms for all perturbations as g = qe 
where hatted quantities are the mode amplitudes, /, m, and k are horizontal and vertical 
wavenumbers respectively, and A is the growth rate. This procedure yields a cubic equation 
for A in terms of the wavevector {l,m,k) and of the non-dimensional parameters: 



A_ 

where K"^ 



= f + m? + k^. 



{WT^) - + ^^') + ^0 '(A + i^^) = , (12) 



It can be shown that the most unstable mode always occurs when the vertical wavenum- 
ber, k, is equal to zero. This corresponds to an "elevator mode", similar to the ones found 
in the related problem of fingering convection and of homogeneous Rayleigh-Benard convec- 
tion. Furthermore, without loss of generality, we may select / or m to be zero by reorienting 
the Cartesian frame so that the fastest growing modes can be described using one horizontal 
wavenumber only. The cubic equation ( !T2|) then becomes: 

A 



Pr 



+ l']iX + niX + Tl') - (A + Tl') + i?o '(A + /^) = . 



(13) 
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There are typically one real and two complex roots of this equation. The real root is always 
negative, achieving a maximum of zero when / = 0. However, the complex root yields a 
positive growth rate when the inverse density ratio lies in the interval Rq^ G (l, f^^)- The 
oscillation frequency is close to the buoyancy frequency. The unstable modes can therefore 
be viewed as overstable gravity waves, as mentioned earlier. 

By maximizing the real part of A over all horizontal wavenumbers, we can find the 
most unstable mode in the system. Figure [T] shows its wavenumber and growth rate for 
Pr = r = 1/3. By and large, the most rapidly growing mode has wavelength of the order 
of 20d regardless of Rq^ for the parameters selected, which implies a lengthscale of a few 
kilometers as described in §2.11 

0.09 
0.08 
0.07 
0.06 
2 0.05 

CD 

Dt= 0.04 
0.03 
0.02 
0.01 


1 1.2 1.4 1.6 1.8 2 2.2 

Fig. 1. — Growth rate (solid line) and horizontal wavenumber / (dotted line) of the most 
rapidly growing gravity wave (i.e. the double- diffusive mode of instability), as a function of 
the inverse density ratio, for Pr = r = 1/3. All quantities plotted are in the units used in 
this paper (see §2.11) . 
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2.3. Ro vs. Rq 

As mentioned previously, it is standard in oceanographic studies of double-diffusive 
convection to use the parameter Rq^ = (3noz/aToz, while Rq = aToz//3noz is used in studies 
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of fingering convectioiJl]. The reason for this change in convention is to emphasize the 
symmetry between the two regimes, which is apparent in Figure [2l Indeed, in that case, one 
can write that overturning convection occurs for 

Rq < 1 fingering regime, 

Rq^ < 1 diffusive regime, (14) 
the double- diffusive regime occurs for 

^ ^ Rq — fingering regime, 
r 

Pr + 1 

1 < i?Q ^ < diffusive regime, (15) 

Pr + r 

and finally that the system is stable for 

Rq > — fingering regime, 

r 

Pr + 1 

> diffusive regime. (16) 

Pr + r 

The symmetry is even more striking at low Prandtl number, where the critical value for 
complete stability is much larger than one in both cases. For the diffusive regime, this implies 
that the system is fully convective up to the Ledoux-stability criterion, and then double- 
diffusively convective between Ledoux-stability and Schwarzchild-stability. A summary of 
the various regimes of instability, for the compositionally homogeneous case, for the fingering 
case and for this diffusive case is presented in Figure O 



3. Numerical Experiments 

3.1. Description of the experiments 

We solve the governing equations and boundary co nditions presented in ^2. II using a high 
performance spectral code developed by S. Stellmach (ITraxler et al.ll2010bt IStellmach et al. 
2010l : iTraxler et al.ll2010al ). specifically designed for the study of fingering convection in the 



oceanographic context. The code can be used "as is" to model double-diffusive convection 
simply by reversing the sign of both background gradients. Since it was developed for 



^For added confusion, studies of double-diffusive convection often call /S^Qz/aTQ^ the density ratio and 
denote it as Rq although we will not use this naming convention here. 
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studying oceanic convection, the code does not consider entropy separately from temperature. 
As a result, it intrinsically assumes that Tg^ = 0, so that Rq = aTo^/ f3fioz- By virtue of 
being non-dimensional, however, the results can nevertheless straightforwardly be applied to 
the astrophysical context simply by interpreting Rq as being defined by flTT]) . 

Our goals are threefold: (a) to characterize transport by homogeneous double- diffusive 
convection (i.e. in the absence of layers), (b) to determine if, under which conditions, and 
through which process thermo-compositional layers may form and (c) to characterize trans- 
port by layered double- diffusive convection if appropriate. 

For this purpose, we ran a sequence of exploratory numerical experiments setting Pr = 
1/3 and r = 1/3. We selected Pr and r below unity to be in the "low Prandtl number, 
low diffusivity ratio" regime, but not too small so that we could remain in a numerically- 
tractable region of parameter space. While it is po ssible to run simulations for lower values 



of these parameters in 3D (see lTraxler et al.ll2010al . for example), these are computationally 



much more demanding. In this first analysis, we wanted to be able to run a uniform set of 
simulations across the whole instability range, and integrate some of them for a significant 
length of time to observe the layer formation and merger process if and when it occurs. We 
discuss the applicability of our results to lower Prandtl number and lower diffusivity ratio 
environments in ^ 

In order to get statistically meaningful measurements of the turbulent fluxes in this 
system, to address points (a) and (c) raised above, we must ensure that the computational 
box contains at least a few wavelengths of the most unstable mode in the horizontal directions 



( iTraxler et al. 1 1201 Obi ). As found in Figure [H these are of the order of about 20d so we select 
a domain size with = Ly = lOOd in all simulations. We use an aspect ratio of 1 and 
choose Lz = 100, which corresponds to Ra-r = 10^. For comparison, we also ran a series 
of simulations with narrower domains {L^ = Ly = 50d) at Ra-r = 10^ for all Rq^, as well 
as one taller-domain simulations using = 178d (equivalently, Rar = 10^) for Rq^ = 1.2. 
See Table 1 for a comparison of these runs. In dimensional terms, the domain height is of 
the order of a few tens of kilometers using the estimates presented in §2.11 for typical stellar 
interiors, which is much smaller than a pressure or density scaleheight and therefore fully 
justifies the use of the Boussinesq approximation. 

In terms of spatial resolution, we use a sufficient number of Fourier modes in all simu- 
lations to resolve the typical size of the composition and velocity perturbations (which are 
roughly of the same size when r = Pr). Similar runs with different spatial resolutions are 
presented in Table 1, and show consistent results. 
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Table 1: Summary of runs performed and measured turbulent fluxes in the homogeneous 
phase. 
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3.0±.4 


4.7±1.3 


.57±.09 


1.35 


96 


96 


420 


900 


1.9±.3 


2.5±.8 


.51±.ll 


1.35^") 


96 


96 


430 


1320 


1.8±.3 


2.2±.9 


.49±.15 


1.5 


96 


96 


650 


1450 


1.5±.2 


1.8±.5 


.51±.09 


1.5(") 


96 


96 


600 


1530 


1.5±.2 


1.7±.8 


.49±.17 


1.6 


96 


96 


800 


1100 


1.4±.l 


1.5±.3 


.52±.07 


1.6(") 


96 


96 


750 


880 


1.4±.l 


1.5±.4 


.52±.09 


1.85 


48 


96 


2000 


2500 


1.17±.03 


1.2±.l 


.57±.04 


1.85^") 


48 


96 


1940 


2500 


1.19±.05 


1.2±.2 


.57±.05 


2.1 


48 


96 





2500 


1 


1 


.6306 



Note. — For each value of the inverse density ratio Rq , the second and third cohimns show the horizontal 
and vertical number of grid-points used. The times tstart and tend define the interval over which the turbulent 
fluxes, measured via Nut and Nu^, were averaged. The quantity 71 is the primary way of measuring the 
total flux ratio 7*°*; see text for detail. Errors quoted indicate the rms fluctuations of NuT(t), Nu^(t) and 
1/71 (t) about the respective means. Unless otherwise specified (see footnotes), simulations were run at 
domain sizes of lOOd x lOOd x lOOd. Note that the simulation for Rq^ = 2.1 was integrated for significant 
time, but by the end of the run the perturbations had not grown to high-enough amplitudes to induce a 
statistically significant heat or compositional flux. 

"50(1 X 50d X lOOd 
X Ud X 17M 

<=89d X 89d X 178d 
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3.2. Qualitative description of the results 

We ran a first set of cubic-domain simulations {Lx = Ly = Lz = lOOd) for Rq^ varying 
across the whole instability range, which for our selected parameters corresponds to 1 < 
Rq^ < 2.167. We found that the initial behavior of the system is qualitatively similar 
for all density ratios: the perturbations first grow exponentially, and then saturate into a 
homogeneous turbulent state. 

The initial exponential growth is well-approximated by linear theory. This is illustrated 
in Figure [3] for example, which shows the early temporal evolution of the rms velocity for the 
Rq^ = 1.2 run, and compares it with the growth of the fastest-growing mode according to 
linear theory. The fit is very good - the small discrepancy at early times can be attributed 
to the fact that more than one mode are excited, but the most rapidly growing mode then 
quickly takes over. We confirmed that the initial instability is independent of the Rayleigh 
number (i.e. the domain height) by comparing these results with those of a taller-domain 
simulation (L^ = 178) at the same Rq^- 

The perturbations saturate once nonlinear effects become important. We find that the 
level of saturation of the turbulence is also independent of the domain height (alternatively, 
of the Rayleigh number), as shown in Figure O However, it depends sensitively on the value 
of the density ratio. Figure H] shows the temporal evolution of the thermal Nusselt number 
Nut as a function of the inverse density ratio. A Nusselt number is the ratio of the total 
flux (diffusive -|- turbulent) to the diffusive flux, so we define the thermal Nusselt number as 

Nut = = l+<wT> , (17) 

where the angular brackets denote a spatial average over the entire domain. We also define 
the equivalent compositional Nusselt number Nu^ as 

Nu^ = = 1 H < wfi > . (18) 

In each case, in the second expression the turbulent fiuxes are expressed in non-dimensional 
form recalling that Tq^ < and /ioz < while T and /i are non-dimensionalized using \Tqz\- 

In all simulations presented in Figure HJ the thermal Nusselt number increases expo- 
nentially until saturation, and remains approximately constant during the early saturated 
phase. After saturation, however, simulations which were run using a lower Rq^ behave in a 
fundamentally different way from those at higher Rq^. In latter case, (for i?o ^ > 1-35), the 
Nusselt number at saturation remains statistically steady for the entire duration of the run. 
By contrast, for Rq^ < 1.35, the Nusselt number later continues to increase. 
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When visualizing the results, we find that this second increase in the turbulent transport 
properties of the system corresponds to the formation of well-mixed fully convective layers 
separated by thin stably stratified interfaces (see Figure O for example). The Nusselt number 
continues to increase as the layers merge, until a single layer is left. We have therefore 
established that layers can indeed form in low-Prandtl number double- diffusive convection, 
and that, as in the high-Prandtl number regime, a layered system transports heat more 
efficiently than a homogeneous system with the same overall temperature and compositional 
gradient. We now study both the homogeneous phase and the layered phase in more detail 
in §l]and ^respectively. 



4. Homogeneous Double-Diffusive Convection 

We focus here on the homogeneous phase, prior to the formation of the first set of layers, 
and measure the transport properties of the turbulence via the respective Nusselt numbers 
defined in (I17p and (1181) . Note that the time period between the initial saturation of the 
double-diffusive instability and the onset of layer formation, when it occurs, varies with R^^ 
(see Figure H]). Table 1 indicates, for each simulation, the time interval over which the system 
is in this homogeneously turbulent phase and during which we average the instantaneous 
Nusselt numbers. 

The mean Nusselt numbers thus extracted are presented in Table 1 and illustrated in 
Figure |6^. The errors quoted denote the rms of the fluctuations around the respective means. 
For most values of Rq^, we ran a series of simulations with different resolution, or different 
box size, or both. The measured Nusselt numbers are always consistent within the errorbars. 

As expected, turbulent mixing is negligible close to marginal stability, i.e. when Rq^ — t- 
(Pr + l)/(Pr + r). It increases as Rq^ decreases through the instability range, and grows 
rapidly close to the onset of overturning convection (i.e. as Rq^ — 1). However, we flnd 
that it remains fairly weak, with Nut of the order of a few and Nu^ of the order of ten, 
even for our lowest Rq^ run {Rq^ = 1.1). By contrast, homogeneous Rayleigh-Benard 
convection in the absence of compositional gradient, in the same parameter regime (Pr 
= 1/3, Ray = 10*^ and as pect ratio one), wo uld have a thermal Nusselt number of the 



order of several thousands (IGaraud et al.l 120101 ). This shows that while turbulent mixing is 



not negligible in this homogeneous double-diffusive regime, it is nevertheless much smaller 
than that induced by standard convection. Presumably, there exists a very narrow range 
of inverse density ratio close to unity, but above it, across which turbulent mixing rapidly 
but continuously increases towards the fully convective value. This will be the subject of a 
subsequent study. 
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Since this quantity is crucial to the theory of layer formation in the fingering regime 
( lRadkdl2003l : iTraxler et al.ll2010al ). we also compute the so-called total flux ratio 7*°*, defined 
as the ratio of the total buoyancy flux due to heat transport to that due to compositional 
transport^]: 

RqNut 



7tot 



T7tot 

ptot 



Nu„ 



(19) 



There are two different ways of computing 7tot from our numerical results: the first and 
preferable method involves calculating 7tot(^) at every timestep|f| and then taking its mean 
during the homogeneous turbulent phase. The second involves using the measured mean 
values of Nuy and Nu^ directly into (IT^ . These two versions of 7tot are denoted as 71 and 
72 respectively, in Figure [6]3. They are consistent within the errorbars. 

Note that we are actually plotting the function 7tot(-Ro^) Figure [6]d. The reason for 
showing it rather than 7tot(-R(7^) will be revealed in §5.21 and is related to the layer formation 
mechanism. As we shall see, the fact that 7^"^ is a decreasing function of R^^ in the interval 
[1, 1.35], i.e. the same interval in which layers are observed to emerge, is not a coincidence. 



5. Layered convection 

As discussed in §3.2[ in all simulations with R^^ < 1.35 we find that the system does 
not remain for long in a state of homogeneous double-diffusive convection, but spontaneously 
develops thermo-compositional layers instead. We study this process in more detail in this 
section, focusing on the Rq^ = 1.2 run. We choose this value of Rq^ rather than one closer 
to unity, based on the results of Figure HI Indeed, for Rq^ = 1.2 layer formation and mergers 
are "slow enough" to be studied, but proceed much more rapidly for lower Rq^. 



5.1. General considerations 

The results presented in Figure H] are for a cubic box with height = 100. In that run, 
we observe that two layers initially form, then merge into one a little while later. In order 



^Note that and are dimensionless here. The interpretation of 7*°* as a buoyancy flux is more 
apparent when we go back to the dimensional quantities and ^^ot.dim^ since = f!^^''^^"^ / kt\Toz\ 

and = F^°*^'i™/(a//3)KT|ro.|. 

■^To be precise, we calculate 7tot(*) ^t every timestep, take the average of this function (shown in Figure 
[5)3), and then take its inverse to get the mean 7tot- The reason for doing this in two steps is that there are 
very occasional events where Nu^(t) < 0, at which point 7tot(^) ~^ 00 but the inverse remains well-defined. 
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to minimize the influence of tlie finite domain height on the layer formation and merger 
process, we use from here on the taller-domain simulation, for which = 178 (Ray = 10^). 
The appearance and successive merger of layers observed in that run is shown in Figure 
O We find that twice as many layers initially form in this nearly-twice-as-tall domain, and 
appear roughly at the same time as they did in the cubic box run. This shows that layer 
formation depends only on local processes, and knows about the thermal scale d rather than 
the domain scale. The initial layer height, in each case, is of the order of 45(1 — 50d. 

Figure [7] shows the evolution of the thermal and compositional Nusselt numbers in 
the tall-domain simulation. We see quite clearly the stepwise increase in transport which 
accompanies the layer formation and successive mergers: there are five fairly well-defined 
plateaus, corresponding to the homogeneous phase, (up to about t = 1000, see Figure Ek), 
the four-layer phase (up to about t = 1200, see Figure Eb), the three-layer phase (up to 
t = 1450, see Figure [St), the two-layer phase (up to t = 1650, see Figure [Sji), and the final 
one-layer phase, see Figure [5^. It is also interesting to note how the heat and compositional 
fluxes follow each other closely throughout the entire simulation. 

A useful way of s tudying the formation and structure of the layers was presented by 



Stellmach et al.l ( 120101 ). and consists in looking at Fourier modes of the density perturbation: 

~p{x, y, z,t) = J2 A,m,fc(t)e^'"+^™^+^'^ , (20) 



Lm,k 



where k is an integer multiple of 2tt/Lz and similarly for / and m. These Fourier modes are 
straightforwardly extracted from our numerical solutions since our code is spectral. 

By definition, the Pco.fc modes are the vertical Fourier modes of the horizontally-averaged 
density profile. Since thermo-compositional staircases are also density staircases (i.e. with a 
nearly uniform density within the layers separated by sharp pycnoclines), a staircase with n 
layers has a dominant vertical wavenumber = 2TTn/Lz. This is seen most clearly in Figure 
[HI which shows the norm of po,o,fc for the four gravest non-zero modes, and illustrates how 
ki, ks, ^2 and ki successively take over as the dominant mode as the layers form and merge. 



A rather striking feature of Figure [HI however, is that the A;4 and k^ modes actually 
begin to grow as early as t = 500, long before the layers appear in visual inspection of 
the temperature and composition field (as in Figure [5] for example). An even more striking 
result is that this growth is well-approximated by an exponential. This strongly suggests 
that layer formation arises through a secondary linear instability of homogeneous double- 
diffusive convection (see next section) rather than throu gh stochastic overturning even ts of 
the growing gravity waves, as is commonly assumed (e.g. IStevensonI Il979t ISpruitlll992l ). 
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We now study in more detail, successively, the layer formation process, and the evolution 
of the fluxes through the staircase as the mergers proceed. 



5.2. The 7-instability 
5.2.1. The 'J— instability of fingering convection 



Recently, significant progress has been made in understanding the spont aneous forma- 
tion of layers in fingering convection, thanks to the groundbreaking work of iRadkd (120031 ) 
in the oceanographic context. Radko discovered that homogeneous fingering convection in 
that regime is linearly unstable to a secondary large-scale instability, which takes the form 
of exponentially-growing horizontally-invariant perturbations in the density (equivalently, 
temperature and composition) profile. The perturbations grow through positive feedback 
between the perturbed stratification and the modulated fingering fluxes. Upon reaching a 
critical amplitude, the modulated density profile becomes unstable to direct overturning con- 
vection, and the system rapidly transitions into a fully-formed, regularly-spaced staircase. 
A sufficient condition for this instability to occur is that the ratio of the total buoyancy 
flux due to heat to the total buoyancy flux due to salt, the quantity referred to as 7tot in 
§11 should be a strictly decreasing function of the density ratio. For this reason, this new 
instability was called the 7— instability and the associated perturbations, the 7— modes. 



Radko's theory was validated first thro ugh two-dimensional s in aulations (IRadkd 120031) . 
and more crucially through 3D simulations (jStellmach et al.ll2010h . IStellmach et al.l (120101 ) 
were the first to find spontaneous layer formation in 3D simulations of fingering convection 
in the "oceanic" parameter regime. They analyzed the growth and structure of the emergent 
layers using the vertical Fourier modes of the density profile, as we have done in the previous 
section. They found, as we do in Figure [HI that the Fourier mode which corresponds to the 
number of layers of the emerging staircase began to grow long before the layers form, and 
that its growth rate could be predicted very accurately by Radko's 7— instability theory (see 
their Figure 6). 



Traxler et al.l (l2010al ) extended Radko's theory for layer formation in fingering convec- 



tion to a parameter regime more relevant of the astrophysical context. Their results suggest 
that at low Prandtl number and low diffusivity ratio 7tot is always an increasing function 
of the density ratio, so that 7— modes are stable. They concluded that spontaneous layer 
formation is unlikely in astrophysical fingering convection. 
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5.2.2. The 'y— instability of double- diffusive convection 



Radko's theory is quite general, and can be applied with only minor modifications to the 
case of double-diffusive convection. Let us consider a system in a homogeneously turbulent 
state, with a background density ratio Rq. We know through the results of §l]that this system 
drives a non-zero total vertical heat flux and a non-zero total vertical compositional flux 



FT with 



TTltot 



Nur(/?o) and F. 



tot 



TTltot 



(21) 



7tot(^o) 

As long as the background is homogeneous, however, and -F^°* are constant in the 
domain, and do not affect the temperature or chemical composition of the system. 

Let us now assume that this homogeneously turbulent system is modulated by a large- 
scale, horizontally-invariant perturbation, so that the horizontally-averaged temperature and 
mean-molecular weight profiles, T and /i can be written as 



T{z,t) = fe 



ikz+M 



(22) 



and similarly for /i. These large-scale perturbations change the local density ratio, which we 
now write as Rp{z, t) = Ro + R'{z, t). Since the turbulent fluxes are functions of Rp, as shown 
in Figure Et, then the perturbations also induce a spatial modulation of the turbulent fluxes. 
In adequate circumstances, the convergence/divergence of the modulated fluxes reinforce the 
original temperature and compositional perturbations, and close the feedback loop. We show 
in Appendix A that the growth rate of the perturbations. A, is the solution of the following 
quadratic: 



A2(l-^)+Nuo(l-Aii?o) 



AiA;^i?oNuo = 



where 



Ai=Rc 



tot J 



dRr 



A2 = Ro 



dNuT 



Ro 



dRn 



Ro 



Nuo = Nut(-Ro), 7o = 7tot(-Ro) • 



(23) 

(24) 
(25) 



This q u adrat ic is exactly the same as that of iTraxler et al.l (j2010al ) and by proxy that of 
Radkd (120031 ) provided his 7 is interp reted as the total flux ratio 7tot, see Appendix A for 
detail. As discussed by iRadkd (120031 ) there is a positive real root when Ai > 0, i.e. when 
7tot decreases with increasing density ratio, or alternatively, when 7^"^ decreases with Rq^. 
Exactly as in the case of flngering convection, a necessary condition for layer formation in 
double- diffusive convection is that '-/tot should be a decreasing function of Rq. 
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Based on the results of Figure |6]3, we can now straightforwardly explain, thanks to this 
theory, the dichotomy between the cases with Rq^ > 1.35, for which we do not expect (and 
do not observe) layer formation, and the cases with Rq^ < 1.35 for which we do. The case 
Rq^ = 1.35 is unclear given our measurement errors on 7tot- Furthermore, since the mode 
growth rate increases with Ai, we can also explain, at least qualitatively, why the staircase 
forms much more rapidly at lower Rq^ (i.e. because the ■jtot curve is steeper in that regime). 
In what follows, we now compare theory and simulations more quantitatively. 



5.3. Comparison of the 7— instability theory with numerical results 

The tall-domain simulation described in §5.11 shows the emergence of a four-layer stair- 
case. Based on the discussion of the 7— instability of §5.2[ we need to compare the growth 
rate of the mode shown in Figure |8]to the solution of ( |23|1 with k = k^. In order to do this, 
we first estimate the coefficients Ai and A2. The derivatives of Nut and 7tot with respect to 
Rp, at Rp = Rq = are calculated using the fluxes obtained in simulations at neighboring 
values of the inverse density ratio. We flnd that 

Ai = 0.453 , A2 = 12.9 . (26) 

We also use the value of Nuq (and associated errorbar) given in Table 1 for R^^ = 1.2, and 
the average of the two 7tot values for 70. 

We find that the dominant 7— mode has a dimensionless growth rate 

A(fc4) = 9.42 X 10"^ . (27) 

Similarly, we find that 

Aiks) = A(A;4)-r| = 5.30 x 10"^ . (28) 

Figure El compares these theoretical predictions with the numerical results for the ^4 and 
fcs modes. We find that they over-predict the growth rate of the dominant k^ mode, but 
provide a fairly accurate estimate of the growth rate of the k^ mode. 

The fact that the ^4 mode grows somewhat slower than predicted is actually expected, 
and leads us to discuss an important caveat of the 7— instability theory. The derivation of 
(!23|) is fundamentally based on an assumption of scale separation between the small-scale 
turbulence, which drives the heat and compositional fiuxes, and the large-scale temperature 
and compositional perturbations T and /2. However, solutions of (!23|) satisfy the similarity 
law A (X k"^, implying that for any given 7— mode there exists, in theory, another more rapidly 
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gr owing one wit h smaller wavelength. This "ultra-violet catastrophe" problem was discussed 
by lRadkd ( 120031 ). and is clearly an artefact of the model. In practice, the 7— instability theory 
should only be applied for k significantly smaller than the wavenumber of the fastest growing 
mode of the primary instability, in order to satisfy the required separation of scales. 

It is therefore interesting to see that the mode which emerges as the dominant 7— mode in 
Figure [8] has a wavelength of about 45(i, which is only about twice as large as the wavelength 
of the most rapidly growing primary gravity wave (about 20d, see §2.2p . As such, it appears 
to be the smallest-scale mode for which the 7— instability theory remains applicable, since 
it still grows exponentially, albeit with a gro wth rate somewhat slow er than predicted. This 
behavior is very similar to the one found by IStellmach et al.l ( 120101 ) in the case of fingering 
convection. It is reassuring to see that, by comparison, the larger-scale mode grows as 
well, and this time with a growth rate which is much closer to the predicted one. 



5.4. Staircase formation 



The dominant 7— mode grows in amplitude until it causes regularly spaced local in- 
versions in the density profile (see Figure |9]). The critical amplitude pcrit, for which the 
total density gradient poz + dp/dz = somewhere in the domain, depends on the mode 



wavenumber and on the background density ratio (e.g. IStellmach et al.l 120101 ): 



Rn 



Pcrit 



(29) 



As the mode grows beyond this amplitude (shown as a horizontal line in Figure [S]), 
progressively larger regions of the domain are unstable to direct overturning convection, and 
a fully-formed staircase rapidly appears. 

Once layers have formed, the 7— instability theory no longer applies and the mode stops 
growing. The subsequent evolution of the staircase through mergers is caused by subtle 
differences in the fluxes through the interfaces and through the convective regions. While a 
full study of the merger dynamics is beyond the scope of this paper, we present in the next 
section an analysis of the heat and compositional fluxes though the staircase as a function 
of the layer height. 
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5.5. Layer mergers and flux increase 

The time of eacli merger rougtily corresponds, in Figure [H to tlie point where the 
amplitude of the kn mode overtakes that of the kn+i mode. Shortly afterward, the new 
dominant mode stops growing, at which point the staircase has reached a new temporary 
n— layered equilibrium. We can then measure the transport properties of the system while 
it is in this n— layered state. The results are presented in Table 2 and illustrated in Figure 

m 

Assuming that the sum of the thicknesses of all the interfaces is small compared with 
the height of the domain, we deduce the layer height HL{n) = L^/n = 178d/n. We can then 
construct a Rayleigh number based on the layer height rather than the domain height, 

Ra. = (^)' . (30) 

Figure [10] shows the variation of Nu^ with Ra^. Rather remarkably, we find that within 
the errorbars the heat flux through the system is more or less consistent with the standard 
scahng law for con vection between two bounded plates separated by a distance (see 



Garaud et al.ll2010l . and references therein). 



Nut ^ OMRa]!^ . (31) 

Note that no fitting was involved here. This scaling would be consistent with an interpreta- 
tion of the interfaces as impermeable, diffusive boundary layers. 

The compositional Nusselt number is about twice the thermal Nusselt number, in the 
layered phase and in the homogeneous phase. A possible explanation for this result is the 
following. The turbulent flux ratio 7 =< wT > / < wjl > is ty pically of order unity for 



turbulence induced by double- diffusive instabilitiecl ( iRadkol 120031 ). As a result, using (fT7|) 
and ( !T8|) we have: 

Nu. ~ 1 + — < >= 1 + — (Nut - 1) = 1 + 2.5 (Nut - 1) (32) 
r r 



"'A plausible reason for this is the following, as argued bv lRadkd (|2003l ). As Rq^ 1, the compositional 
field acts more and more like a passive tracer. As a result, the turbulent diffusivities for heat and composition 
tend to one-another, and since i?o = 1, so do the induced turbulent buoyancy fluxes. As R^^ — > (Pr+r) / (Pr+ 
1), 7 tends to one for a different reason. Since the instability is driven by the conversion of potential energy 
into kinetic energy, near the marginal stability limit this available potential energy must vanish. Hence, the 
turbulent buoyancy flux due to heat must exactly equal that due to composition for the potential energy 
available to vanish. Finally, since 7 has to tend to one in both limits, it cannot deviate significantly away 
from one inbetween. 
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for our selected parameters. For large enough Nusselt numbers, this implies Nu^ ^ 2.5Nut. 
Figure [10] shows f p2|) in comparison with the data. Again, the fit is satisfactory within the 
errorbars, altho ugh not quite as comp elling as that of Nu^. Note that a similar argument 
was invoked by iTraxler et al. J2010ah to explain the relationship between their measured 



scaling laws for Nut and Nu^ in fingering convection. 



It is also interesting to compare our numerical results with the work of ISpruitI ( 1l992l ). 



who proposed the following parametrization for heat transport by layered convection: 

Nur ^ 0.5(PrRai)^/^ . (33) 

His estimate of the turbulent compositional diffusivity (see his equation (44)) yields the 
following compositional Nusselt number: 

Nu^ oc t-^/^^Nut = Rot-^^^Nut (34) 

in this Boussinesq case, where the proportionality constant is of order unity. Both esti- 
mates are shown in Figure [10] as well, and are also more-or-less consistent with our flux 
measurements within the errorbars. This time, Nu^ seems to be better accounted for than 
Nut. 

Further simulations will be needed to determine which (if any) of these scaling laws best 
explains the data. Wider domains will be needed to improve our statistics and reduce the 
variability of the mean fluxes. In addition, taller domains will be necessary to allow for larger 
layer heights. But beyond the details of the power laws or the prefactors, our simulations 
clearly show that the overall transport through a staircase depends only on the layer height, 
and that the heat and compositional fluxes are roughly proportional to each other for given 
Rq^, Pr and r. 



6. Discussion and conclusion 

6.1. Summary of the results 

In this work, we have studied a set of numerical simulations of double-diffusive convec- 
tion, in a triply-periodic domain, for Prandtl number Pr = u/kt = 1/3 and diffusivity ratio 
r = Kh/kt = 1/3. We have explored the entire instability range, varying the inverse density 
ratio Rq^ between 1 (the onset of direct overturning convection) and (Pr + l)/(Pr + r) (the 
marginal stability limit). Our simulations were performed in a "small" domain spanning, 
in the horizontal direction, about five wavelengths of the fastest-growing double-diffusive 
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mode (i.e. = Ly = lOOd where ci is a thermal diffusion lengthscale), and in the vertical 
direction, = lOOd or = I78d depending on the runs. 

In all cases we initialized a double-diffusively unstable system with infinitesimal per- 
turbations, and found that these first grow exponentially according to linear theory, then 
saturate into a state of homogeneous double-diffusive convection. In that state, the turbulent 
contribution to thermal and compositional transport is significant but much smaller than 
that expected from standard convection, ranging from 5-10 times the diffusive rate near the 
onset of direct convective instability, and rapidly dropping towards zero as Rq 
towards marginal stability (see [0]). 



mcreases 



For small Rq^, however, the system does not remain in this homogeneously convecting 
state. Instead, thermo-compositional layers rapidly appear, and transport through the sys- 
tem strongly increase s. We showed that the layer format i on process is governe d by Radko's 
7— instability theory (IRadkol 120031 : IStellmach et al.ll2010l : iTraxler et al.ll2010al ). both quali- 



tatively and quantitatively. In particular, it explains why our simulations with Rq^ < 1.35 
transition into layers while those with Rq^ > 1.35 do not. The key factor is the variation of 
the total buoyancy flux ratio 7tot with density ratio (see §5.2p : layers can only form when 
7tot decreases with _Ro- 

In the layered phase, we found that the flux through the staircase depends sensitively 
on the mean layer height Hi. Given the large variability of the measured fluxes during the 
layered phase, our results are roughly consistent both with Spruit's theory (ISpruitlll992l ) and 
with heat transport between two solid plates (as given by equation (13T]) ). Further simulations 
will be needed to help distinguish between these two possibilities - or perhaps suggest an 
alternative one. Finally, note that in our small-domain simulations, the mergers always 
proceed until a single layer is left. In that sense, the dynamical evolution of the system is 
always eventually influenced by the domain size. 



6.2. Discussion of the applicability of our results to real systems 

Our initial goals were threefold: (a) to characterize transport by homogeneous double- 
diffusive convection (i.e. in the absence of layers), (b) to determine if, under which conditions, 
and through which process thermo-compositional layers may form and (c) to characterize 
transport by layered double-diffusive convection when appropriate. 

To answer part (a) in detail, a much larger number of simulations will be needed, using 
progressively smaller Prandtl numbers and diffusivity ratios. These are the subject of an 
ongoing investigation. We hope to find similar scaling laws for Nu^ and Nu^ as functions of 
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Pr and r as the one found by iTraxler et al.l (l2010bl ) for fingering convection 



By contrast with the case of fingering convection, however, we now know that thermo- 
compositional staircases can form spontaneously from double- diffusive convection. Radko's 
criterion for layer formation, namely that 7tot should decrease with Ro, is the answer to 
part (b) of our goals, but does require knowledge of the function 7tot(-Rp) to be applied in 
practice. The latter must be determined separately for each parameter set (Pr, r). 

Finally, our findings have provided some insight into part (c). Since transport through 
a staircase depends on the layer height only (for given fluid parameters and overall stratifi- 
cation), the problem shifts to estimating actual layer heights in astrophysical objects. The 
layers we observe in our simulations have a strong tendency to merge, which suggests two 
possible outcomes: the mergers continue indefinitely, until the scale of the equilibrium layers 
is commensurate with the system size; or the mergers stop for other reasons (see below), with 
an equilibrium layer height significantly smaller than the system size. It is of course crucial 
to know which of these two scenarios is correct, as they imply vastly different transport rates 
through the staircase. 

Unfortunately, our simulations were not able to provide a definitive answer to this 
question. In the two cases studied, the final layer height was equal to the domain height, but 
this should not be interpreted as a result in favor of the first scenario since this could simply 
mean that our domain was too small to "contain" the intrinsic equilibrium layer height of 
the second scenario. 



Radkd (120051 ) proposed a theory supporting the idea that mergers stop before layers 
reach the system size, and deduced a means of estimating the equilibrium layer height. 
Starting from an initial staircase with uniform "jumps" in temperature and chemical com- 
position across the interfaces, he studied how the staircase evolves if it is perturbed slightly, 
by making some of the jumps larger and some of the jumps weaker. He concluded that the 
staircase is unstable to mergers if the total buoyancy flux ratio through the interfaces is a 
decreasing function of the density ratio across the interfaces - a criterion very similar to the 
7— instability criterion. 

We have tried to test Radko's merger theory against our simulations, but this has 
unfortunately proven to be difficult. The statistical fluctuations in the measured turbulent 
fluxes (see Figure [7]for example) are too large to detect a significant variation of the interfacial 
fiux ratio as the mergers proceed. We would need a much larger domain to improve the signal 
to "noise" ratio to a point where our results could be compared with his theory. We would 
also need a much taller domain (at least a few times taller than the equilibrium layer height) 
to see if the merger process indeed stops as expected. Finally, we would need to integrate 
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the simulation long enough to establish convincingly that the mergers have indeed stopped. 
Unfortunately, running equivalent simulations in a much larger domain, and for long enough 
to observe the layer formation and merger process, is impossible within current numerical 
limitations. 



6.3. Future prospects 

The preliminary findings presented in this paper still enable us to lay out a clear path 
towards obtaining better parametrizations of mixing by double-diffusive convection in the 
near future, using currently available computational resources: 

• Firstly, we must gain a better understanding of the instability saturation mechanism 
at low Prandtl number and low diffusivity ratio, in order to determine the flux laws 
NuT(Pr, r, Rq) and Nu^(Pr, r, Rq) for homogeneous double-diffusive convection. These 
flux laws are needed to determine when layers are expected to form, and can be used 
"as is" to parametrize mixing otherwise. They can be measured using "small-domain" 
simulations similar to the ones we have presented here, at least for values of Pr and 
T as low as about 0.01 or so. Semi-analytical weakly-nonlinear models will then be 
helpful to guide extrapolations to the much lower parameter values appropriate of the 
astrophysical regime. 

• Secondly, we must gain a general understanding of mixing in the layered case, at low 
Prandtl number and low diffusivity ratio. In order to do this, we need to determine 
how interfacial transport depends on the fluid parameters (Pr, r) and on the interfacial 
density ratio (ie. a density ratio based on the difference in temperature and composition 
across the layers). We must also understand how transport scales within the convective 
layers, as a function of the same parameters but also as a function of the layer height. 
This can be done today using simulations in which a single layer is pre-seeded, to 
bypass the rather lengthy layer formation and merger phases. Using this information, 
we will be able to te st the basic fl ux laws which are central to Radko's merger theory 



more quantitatively (lRadkdl2005l ). If this theory holds, then one can straightforwardly 
deduce the equilibrium layer height for a given parameter set, and ultimately quantify 
the staircase transport properties. 
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A. Derivation of the 7— instability in the diffusive case. 



Following iRadkd (120031 ) and iTraxler et al.l ( l2010al ). we begin with the general non- 
dimensional governing equations f lTOj) . and average them over several wavelengths of the 
fastest growing mode of the primary instability. We get 
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where R is the Reynolds stress, F^ and Fj^ are the total heat and compositional fluxes 
respectively, and T, /i, p, and u now denote large-scale fields only. 



The 7— instability drives horizontally-invariant perturbation with zero mean flow ( iRadko 



20031 ). We can therefore neglect the momentum equation, set u = 0, and ignore all horizontal 
derivatives. The mean temperature and composition equations simplify to: 
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Finally, we assume that Nuy, and 7tot depend only on the local value of the density 
ratio Rp. Note that Rp is no longer constant, but varies with 2; as a result of the large-scale 
background temperature and compositional perturbations, as 

a (To, + Tf"^) 



Rn 



Ro(l-T, 



(A5) 



' /3(;Uo. + /if l-Rol^z ' 

where, for clarity, we first expressed Rp as the ratio of dimensional quantities and then as 
the ratio of non-dimensional quantities. 

We now linearize equations (lA2p and (1A3P around a state of homogeneous turbulent 
convection in which T = + T',fj, = + fj,', and Rp = Rq + R' where linearization of (lASD 
yields: 

R! = i?o(l -T, + Rof^z) . (A6) 
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Using the fact that F^°^ = Nur(-Rp)(l - T,) and F*°* = -Fy V7tot the hnearized temperature 
equation becomes 



dT' 

~dt 
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while the hnearized composition equation is 
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where we have used the foUowing notation for simphcity: 
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This quadratic is exactly the same as the one obtained in the fingering case. In hindsight, 
this result is trivial, and could be obtained immediately had we allowed ourselves to non- 
dimensionalize T and /i using negative dimensions (in which case the governing equations 
and all definitions are exactly the same as in the fingering case). 
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Fig. 2. — Summary of the various regimes of convective instability. The top hue corresponds 
to the chemically homogeneous case, with = 0. The system is unstable to direct overturn- 
ing convection if V > Vad or equivalently Tqz < Tq^ < 0. The Schwarzchild criterion (dotted 
line) appropriately marks the stability boundary (solid line). In the presence of an unstable 
mean molecular weight gradient (middle line), the region unstable to convective overturning 
extends into the subadiabatic regime, and is stabilized only when V — Vad = (Ledoux 
criterion, vertical dashed line). The Schwarzchild criterion in this regime is not relevant. 
Beyond the Ledoux limit, the system can stillhe unstable, this time to fingering convection. 
In the case where the system has a stable mean molecular weight gradient (bottom line), 
the region of parameter space unstable to overturning convection shrinks according to the 
Ledoux criterion. The system can still be unstable to double-diffusive convection in a subset 
of the interval between Ledoux-stability and Schwarzchild stability. 
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Fig. 3. — Temporal evolution of the non-dimensional rms velocity in simulations with Rq = 
1.2, for two different computational domain heights: Lz = 100 corresponding to Ra^ = 10^ 
(dashed line), and = 178 corresponding to Ra^ = 10^ (solid line). The straight dotted 
line shows an estimate of the early exponential growth based on the growth rate of the 
most rapidly growing mode only (see §2.2p . The two simulations saturate at the same level, 
confirming that the dynamics of the system in the saturated phase are independent of the 
domain size (for large enough domains). 

Table 2: Variation of the Nusselt numbers as a function of the number of layers, for the 
i?o ^ = 1-2, Lz = I78d run. 



n 


Hl 


^start 


'^end 


Nur 


Nu^ 


4 


44.5 


1020 


1240 


10.5 ± 1.7 


20.6 ± 4.1 


3 


59.3 


1240 


1480 


16.4± 2.9 


33.5 ± 8.4 


2 


88.9 


1480 


1660 


25.8 ± 5.0 


55.0 ± 12.6 


1 


177.8 


1660 


2190 


46.0 ± 14.5 


99.0 ± 33.3 



Note. — The first column shows the number of layers, and the second column an estimate of the layer 
depth in units of d, = Lz/n. The times istart and tend mark the interval of the time over which the 
turbulent fluxes, measured via Nut and Nu^^, were averaged. Errorbars represent the rms fluctuations about 
the mean. 
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Fig. 4. — Evolution of the thermal Nusselt number for seven of the simulations presented 
in Table 1. In all cases, Pr = r = 1/3, Ra^ = 10® [L^ = 100), and the aspect ratio is one. 
Nut — 1 is shown to emphasize the exponential growth phase. The results are also staggered 
for clarity, so each curve actually shows /(Nu-r — 1) where the multiplicative factor / is 1, 
10, 100, 1000, 10^ 10^ and 10^ respectively for R^^ = 1.85, 1.6, 1.5, 1.35, 1.2, 1.15 and 1.1. A 
straight horizontal line of the same color in each case marks the point at which Nur — 1 = 1 
for reference (i.e. when turbulent and diffusive fluxes are equal to one another). Note how 
runs with Rq^ > 1.35 remain in a quasi-steady saturated state, while runs with i?o ^ < 1.35 
show a subsequent increase in transport. The Rq^ = 1.35 run was actually integrated until 
t = 2500, but was found to remain at the same saturated level. 
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5. — Volume- rendered visualization of the mean molecular weight perturbation, for 
= 1.2, using the tall-domain simulation (L^ = 178c?). Shown are five snapshots taken 
at different times, (a) in the homogeneous phase at t = 400, (b) in the four-layer phase 
at t = 1100, (c) three-layer phase at t = 1350, (d) two-layer phase at t = 1550 and (e) 
single-layer phase at t = 1850. The color scale is adjusted in each panel to emphasize the 
perturbations, so that /t G [— 0.1,0. Ij^Q^L^ in (a), ft G [—0.25, 0.25] /io^i^z in (b) and (c), 
ft G [-0.4,0.4]fiozLz in (d) and /i G [-0.5,0.5]fiozLz in (e). 
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Fig. 6. — Left: Mean Nusselt numbers in the homogeneously turbulent phase (see Table 2), as 
a function of the inverse density ratio. The Nur measurements are shown as (+) symbols, and 
the Nu^ measurements as (x) symbols. The turbulent contribution to transport, represented 
by Nu — 1 in both cases, goes from a few to zero over the instability range. The Rq^ = 2.1 
point was left out since we measured Nut = Nu^ = 1 in that run. Right: Inverse of the total 
buoyancy flux ratio 7tot(-^o^) measured in our simulations, using two different methods 
(see main text for detail). Note how 7t~J shows a pronounced minimum around Rq^ = 1.4. 
In all cases, the errorbars represent rms fluctuations of the respective functions around the 
mean. 
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Fig. 7. — Evolution of Nu^ (solid line) and Nu^ (dashed line) for the Rq^ = 1.2 tall domain 
run {Lz = 178). This plot shows the step- wise increase in transport through the various 
phases. In the layered phase in particular, the heat transport through the staircase depends 
on the layer height. The horizontal lines indicate the mean compositional Nusselt number 
measured in the four-, three-, two- and one-layer state in the tall-domain run (see §5.51 and 
Table 2 for detail). The solid part of each line indicates the interval of time over which the 
averages were measured. Equivalent lines for Nur are left out to avoid cluttering the plot. 
Note how closely the two curves follow each other throughout the entire run. 
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Fig. 8. — Evolution of the norm of the Fourier coefficient of each of the four gravest Fourier 
modes of the vertical density perturbation profile, |po,o,fcP = Po,o,kPo o k- The figure illustrates 
the initial exponential growth of the ^4 and ks modes, and compares them with the prediction 
from the 7— instability theory (see §5.2p . shown as the same-color straight lines. The mode 
grows until it reaches the critical amplitude for overturning (horizontal black line), see main 
text for detail. Shortly afterward, four equally spaced layers appear (around t = 1000). The 
layers later merge successively, which can be seen here easily as the ks, k2 and ki mode 
respectively take over as the dominant mode. 
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Fig. 9. — Evolution of the density profile of a ^4 7— mode. For simplicity, the density is 
normalized to the density difference across the domain height, and the vertical coordinate 
is in units of the domain height. The dotted line shows the background density profile 
with constant gradient poz, the dashed line the background + perturbation for intermediate 
perturbation amplitude, while the sohd line shows the background + perturbation at the 
critical amplitude for the onset of overturning convection given by (129|) . Note the existence 
of four specific points where the total density gradient poz + dp/dz is zero, which will become 
the middle of the emergent layers. 
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Fig. 10. — Variation of the two Nusselt numbers in the layered phase, for the tall-domain 
Rq^ = 1.2 run, with the Rayleigh number defined with layer height, Ra^. The symbols 
correspond to the measurements presented in Table 2, with Nut shown as (+) symbols and 
Nu^ shown as ( X ) symbols. The errorbars show the standard deviation of the measured fluxes 
about the mean. The various lines cor respond to dif ferent possible scaling laws discussed in 
the main text: the blue lines are from ISpruitI ( 119921 ). for Nur (solid line) and Nu^ (dashed 
line), and the green lines are from fl3Tl) for Nuy (solid line) and fl32|) for Nu^ (dashed line). 



